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We discuss the cosmological evolution of the inflationary gravitational wave background (IGWB) 
in the Randall-Sundrum single-brane model. In braneworld cosmology, in which three-dimensional 
space-like hypersurface that we live in is embedded in five-dimensional anti de Sitter (AdSs) space- 
time, the evolution of gravitational wave (GW) modes is affected by the non-standard expansion 
of the universe and the excitation of the Kaluza-Klein modes (KK- modes). These are significant 
in the high-energy regime of the universe. We numerically evaluate these two effects by solving 
the evolution equation for GWs propagating through the AdSs spacetime. Using a plausible ini- 
\^*; , tial condition from inflation, we find that the excitation of KK-modes can be characterized by a 

f"^ ■ simple scaling relation above the critical frequency / cr i t determined from the length scale of the 

fifth dimension I. The remarkable point is that this relation generally holds as long as the matter 
£SJ , content of the universe is described by the perfect fluid with the equation of state (EOS) p = wp for 

< w < 1. The resultant scaling relation is translated into the energy spectrum of the IGWB as 
Q • Ogw oc f( 3w - 1 )/( 3w + 2 ) f or f > J crit . This indicates that in the radiation dominant case (w = 1/3), 

' the two high-energy effects accidentally compensate each other and the spectrum becomes almost 

the same as the one predicted in the four-dimensional theory, i.e., Hgw oc f°. 

■ 

T— I 1 PACS numbers: 04.30.-w, 04.30.Nk, 04.50.+h, 98.80.-k 

> ■ 

m ■ 

O ! I. INTRODUCTION 

Gravitational waves (GWs) are ultimate probes of the untouched region of the universe. Currently, large scale 
ground-based interferometers (TAMA3OO0, LIGO@, VIRGO Q, GEO6OO0) are enthusiastically trying to detect 
the signals emitted from stellar objects with relativistic motion - supernovae explosions, coalescence of neutron star 
binaries and so on. Among numerous types of GWs, the gravitational wave background (GWB) may possess much 
interesting information on the cosmology, though its detection is expected to be so challenging [5(. Especially, the 
inflationary GWB (IGWB), generated during the inflationary epoch by the quantum fluctuations of the spacetime [|| 
El j| H I, IToL Ull IT^ . IT3I is thought to be one of the most fundamental predictions of the inflationary scenario 
^oHlo. Il7l llSf. Since the history of the cosmological expansion is imprinted in the power spectrum of the IGWB, it 
helps us to understand the extreme ly e arly universe if we can detect the signals by the future space-based experiments, 
• Ch ! such as DECIGO[I3 and BBOH^I^. 

As an ultimate cosmological tool, the IGWB may also be useful to probe the presence of extra-dimensional spaces. 
Recent developments in particle physics suggest that we live in a higher dimensional spacetime. In particular, 
braneworld scenarios have recently attracted much attention theoretically and observationally (for a review, see 
|22|). According to them, we live in a three-dimensional hypersurface (brane) embedded in the higher-dimensional 
spacetime (bulk). While gravity can propagate in the bulk with the curvature scale I, the standard model particles 
are confined to the three-dimensional brane. In the low-energy regime of the universe {HI <gC 1), four-dimensional 
general relativity is successfully recovered and the extra-dimensional effects should be fairly small. On the other hand, 
in the high-energy regime {Hi ^> 1), the localization of gravity is not always guaranteed and a significant deviation 
of the time evolution of GWs from the standard four-dimensional theory is expected. If this scenario is true, the 
spectrum of the IGWB may be significantly modified by the high-energy effects, which can provide a direct probe of 
the extra-dimensions. 

The goal of this paper is to investigate the high-energy effects on the evolution of the IGWB and quantify the 
power spectrum. During the inflationary epoch, the wavelength of the GWs exceeds the Hubble horizon scale due to 
the exponential expansion of the universe and the amplitude of GWs becomes frozen. After the end of inflation, the 
universe enters the decelerated expansion phase and wavelengths of GWs soon become shorter than the Hubble scale. 
When the Hubble horizon scale becomes comparable or smaller than the characteristic size of the extra-dimension, 
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the high-energy effects may significantly affect the evolution of GWs. There are two main high-energy effects : i) 
peculiar cosmological expansion due to the high-energy correction of the Friedmann equation, which enhances the 
spectrum in the high frequency region and ii) excitation of Kaluza-Klein modes (KK- modes) freely escaping from the 
brane to the bulk spacetime, which may suppress the amplitude of the GWs on the brane. While the former effect is 
simply estimated from the expansion rate of the universe, the amount of the latter effect requires the knowledge of the 
wave propagation in the bulk. Focusing on the Randall-Sundrum single-brane model, in which the three-dimensional 
hypersurface (brane) is embedded in the five-dimensional anti-de Sitter (AdSs) spacetime|23j. many authors have 
tried to estimate these effects in an analytic way. However, these analyses were restricted to the idealistic situations 
|2H I2H l2l| or low-energy cases in the Friedmann universe [53, H^, since the analytic study of wave equation is 
generally intractable due to the complicated form of equation as well as boundary condition. Thus, in this paper, we 
numerically solve the wave equation and try to estimate the observed IGWB spectrum. 

Concerning numerical studies, several authors have used different numerical techniques and coordinate systems to 
solve the wave equation of GWs [U HJ HJ IH S H . In our previous studies |29l |30j , numerical simulations were 
carried out in the two types of coordinate systems. One is the Gaussian normal coordinate system in which we found 
the suppression of the amplitude of the IGWB on the brane. Unfortunately, the coordinate singularity appears in 
the bulk and this restricts our analyses to a relatively low-energy scale |29j . On the other hand, another coordinate 
system we used is the Poincare coordinates in which we observed that the two high-energy effects compensate each 
other and the spectrum became same one as predicted in the four-dimensional theory l30|. 

In this paper, we first show that the spectrum of the IGWB estimated in the paper [30| is robust against several 
numerical artifacts, such as the dependence of the regulator brane and initial time. Next, we study the dependence 
of the spectra on the equation of state (EOS) of the universe to clarify how the two high-energy effects change with 
the expansion rate of the universe. 

This paper consists of seven sections. In Sec. 2, we discuss how the spectrum of the IGWB is affected by the 
presence of the extra-dimensions. In Sec. 3, we introduce the Poincare coordinate system and derive the evolution 
equation of GWs. After briefly discussing the details of the numerical simulations and initial conditions in Sec. 4, 
we check our numerical scheme by solving a simple case in Sec. 5. In Sec. 6, we present the numerical results in 
the case that the IGWB re-enters the Hubble horizon during the radiation dominated (RD) universe. Also in Sec. 6, 
the results in cases with other equation of state (EOS) are shown. Finally, Sec. 7 is devoted to the summary and 
conclusions. 

In this paper, we use units in which c = h = 1. M p \ represents the four-dimensional Planck mass/energy. 

II. GRAVITATIONAL WAVE BACKGROUND FROM INFLATION 
A. Standard four-dimensional prediction 

Standard inflation model predicts that gravitational waves (GWs) are generated by the quantum fluctuations of 
spacetime. In this section, we briefly review how to evaluate the power spectrum of IGWB (see also Ref.|T^1. 

In the left panel of Fig 2] a sketch of the evolution histories of GWs with various wavelengths is shown. The vertical 
and the horizontal axes represent the wavelength of GWs and the cosmic time, respectively. In this figure, the solid 
line represents the Hubble horizon scale H -1 , and we have labeled three regimes as "Inflation", "RD" and "MD" 
for the inflationary epoch, the radiation-dominated epoch and the matter-dominated epoch, respectively. During the 
inflation, the universe experiences accelerated expansion and the wavelength of GWs eventually exceeds the Hubble 
horizon scale. Then the oscillatory behavior ceases to exist and the amplitudes of GWs become frozen. After inflation, 
these GWs re-enter the horizon in the decelerated expansion phase (regions "RD" and "MD" ) . Inside the horizon, the 
wavelengths are redshifted and the amplitudes are reduced by the cosmological expansion. Since the horizon re-entry 
time depends on the comoving wave number for each GW mode, the resultant energy spectrum of the IGWB observed 
at present simply reflects the expansion rate at the horizon re-entry time. 

To evaluate the spectrum of IGWB, we first consider the characteristic frequencies of GWB associated with the 
cosmic history. We define three characteristic frequencies according to the standard four-dimensional cosmology : i) 
the lowest frequency fh] ii) the frequency of GWs re-entering the horizon just at the matter-radiation equality time, 
/ oq ; and hi) the cut-off frequency by the inflation, /j n f. First, the largest wavelength of IGWB observed today is 
definitely the horizon length which corresponds to the frequency 

A«M»io-"tt( ratm %. Mpe ). ID 

where Hq denotes the present value of the Hubble parameter. Second, the frequency of large-scale GWs which came 
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into the Hubble horizon at the matter-radiation equality time t eq can be calculated as 

/., = ^.,^.1X10-"H, ( 72k J:. Mpc )(^)' /2 , P) 

where a denotes the scale factor and z the redshift. The subscripts '0' and 'eq' represent the quantities evaluated at 
the present time to and at the matter-radiation equality time t eq , respectively. Finally, the highest frequency observed 
today is determined from the Hubble horizon at the end of the inflation, which can be calculated in the same way as 

i«in f a oqw „ / g inf y/y h y/yi + z^- 1 > 



/inf ~ 27ra oq a ^ ~ L1 GHZ U X 10- 5 M P J 72 km/s • Mpc ) { 3200 J ' (3) 

where Hi n f means the energy scale of the inflation, which is constrained by the COBE observation as iJ in f < 6x 10 _5 M p i 
[l4|. As a consequence, the GWs with /h < / < / C q re-enter the Hubble horizon during the MD phase, while for 
/ eq < / < /inf, GWs re-enter during the RD phase. These characteristic frequencies are shown in Fig. ^ 

Let us focus on the shape of the IGWB spectrum. Conventionally, the power spectrum of GWB is characterized 
by the energy density instead of its amplitude. We introduce the quantity ficw defined as 0] 

M/)^-^, (4) 
Pc a log/ 

where p c = 3Hq/8ttG — 9.8 x 10~ 30 g/cm 3 means the critical density of the universe and pow the energy density of 
GWB. Denoting the characteristic amplitude by h, the above quantity is related to 

^ GW = ( 1.263 x 10-" ) (l4) ■ (5) 

The frequency dependence of the power spectrum can be derived from the fact that the amplitude of the GWs 
evolves as h oc 1/a inside the horizon. Assuming the scale factor evolves as a oc t n near the horizon re-entry time £*, 
the GW amplitude observed today is related to t r as 

/io = -^ocCr T/2 . (6) 

do 

Here /i* denotes the amplitude of GWs evaluated at the time £*. The amplitude /i* is primarily determined by the 
quantum fluctuations generated during the inflation, whose spectral dependence is given by hi oc /™ T (e.g., Sec. 6.5 
of Ref. |3!| ). In a single-field model of slow-roll inflation, the spectral index tit can be expressed by the slow-roll 
parameter e as jit ~ — 2e 36] 35]. In the pure de Sitter expansion, jit = 0. 
In the power-law expansion, the Hubble parameter at t = t* scales as 

H^t- 1 . (7) 

Hence the observed frequency / is related to i* as 

* = ^H 1(xt n- lj (8) 

where k = a»_ff„ denotes the comoving wave number of the GW concerned. From 10 and (J8J, the power spectrum of 
the IGWB becomes 

n Gff «^/ 2 «/^ +,lT . (9) 
Particularly in cases with the matter content in the universe satisfying the equation of state (EOS), 

P = wp, (10) 
the power-law index of the scale factor, n, is rewritten with 
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Then the energy spectrum @ becomes 



(12) 



Now, simply assuming tit = and applying this formula to the standard history of the universe [RD (w = 1/3) phase 
and MD (w — 0) phase] , the energy spectrum of the IGWB observed at present becomes ^3| 



f2 G w oc 




(/cq < / < /inf), 
(/h < / < /cq) 



(13) 



The resultant spectrum in the four-dimensional cosmology is shown in long-dashed lines in Fig. [21 The normalization 
of the spectrum is weakly constrained from the CMB observation by COBE, which leads to f^cw < 10~ 14 |l4|. This 
figure shows that the IGWB widely exists over the frequencies which spans approximately 30 orders of magnitude. 
Furthermore most of the frequency region comes from the GWs which re-enters the horizon during the RD phase. 




cosmic time 



cosmic time 



FIG. 1: Schematic diagrams of the evolution histories of GWs. Left : four- dimensional case, Right : five-dimensional case. 
In the high-energy RD regime of the latter case (denoted by "H"), the p 2 -term in the Friedmann equation changes the time 
dependence of the Hubble parameter. This fact yields a new characteristic frequency / cr it- 



B. High-energy effects in the braneworld cosmology 



In a braneworld scenario, the propagation of gravity may be modified by the presence of extra-dimensional spaces, 
which affects the observed spectrum of IGWB. In this paper, we consider the Randall- Sundrum single-brane model 
|2.'"'J| . In this model, a three-dimensional brane is embedded in five dimensional anti-de Sitter spacetime (AdSs bulk) 
with the curvature scale I. Here we consider the flat Friedmann-Robertson- Walker (FRW) universe on the brane. In 
general, the RS model may possess a black- hole in the bulk, but we do not consider it. 



1. Cosmological Expansion on the Brane 



In the RS model, the Z2 symmetry (mirror symmetry) is imposed on the brane, whose physical meaning comes from 
the S 1 /Z2 orbifold compactification in heterotic M-theory. This symmetry provides the junction condition which gives 
a relation between the extrinsic curvature at the brane and the energy density of matter on the brane. From the 
junction condition, the Friedmann equation on the brane and the conservation law become |37l f38l 13^. l40t l4lL |4^| : 



(14) 
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where k 2 represents the gravity constant defined as k 2 = 8ttG = Sir/M^ and A denotes the tension of the brane. The 
tension is related with the bulk curvature scale £ as k 2 X£ 2 = 6 if we require that the cosmological constant on the 
brane A4 vanishes. For later convenience, we define a dimcnsionlcss energy density normalized by the tension e(t) by 
e(t) = p(t)/X. 

From (|14(l . we can separate the regime into two parts : the high-energy regime (H£ > 1) and the low-energy regime 
(H£ < 1). The critical energy density satisfying the relation HI = 1 is given by 



= V2- 1. 



(15) 



Thus, when e(t) ^> e cr i t , the high-energy correction characterized by the p 2 -term in l|14fl significantly modifies the 
cosmological expansion. 

In the cosmology with perfect fluid, the exact solutions for the scale factor a(t) and the normalized energy density 
e(t) are known. These solutions are expressed as (e.g., j3^, 39]): 



a(t) = (f^*) 



1/3(1+™) 



e(t) 



{3(i+w)t/e+iy 



1 



(16) 



where the scale factor is normalized to unity at the critical energy, e cr i t . In cases with w = (MD), w — 1/3 (RD) 
and w = 1, the energy densities become respectively 



e(t) 



2£ 2 



m 2 + m 

i 2 

8t 2 + AU 
( 2 

1 18i 2 + 6t£ 



for w = 0, 
for w = 1/3, 
for w = 1. 



(17) 



Using the above relations, we will consider how the high-energy regime of the universe can modify the spectrum of 
the IGWB. 



2. High-energy effects on GWs 

As we mentioned in Sec. [I] there are two important effects on the spectrum in the high-energy regime. Let us first 
consider the non-standard cosmological expansion in presence of p 2 -term. From <|16[) . the power-law index of the scale 
factor in the high-energy regime ('IT in the right panel of Fig. ^) is related to the EOS parameter w as 

for t<£. (18) 



3(l+u;) 



On the other hand, at the late-time phase (e <ti e C rit), the power-law index just coincides with the four-dimensional 
result ltTT|l . As a result, the energy spectrum © is modified to 

o GW =( / !!!! fOT/<</ -' (19) 

1^/— for />/ crit , 

where / cr ;t denotes the critical frequency given by 

1 ^crit ^eq 



,/crit — 



2tt£ a cq ao 
5.6 x 10" 5 Hz 



-1/2 / „ \ I/ 2 /I I ~ \ -1/4 ( 20 ) 



0.1mm/ \72 km/s • Mpc / V 3200 



That is, the wavelength at the horizon re-entry time just coincides with the curvature scale, namely, H* = £~ 1 (cf. 
|43|]). Note that the curvature radius £ is constrained to £ < 0.1 mm by table-top experiments on the Newton force 
law 0, . For the frequency / > / C1 it , GWs re-enter the horizon during the RD phase of the p 2 -term dominated 
epoch (right panel of Fig. In the high-energy RD phase, the spectrum of the IGWB is modified to 

fiGW « / 4/3 (/crit < / < /mf), (21) 
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which is shown in short-dashed lines in Fig [21 Additionally, the inflationary cut-off frequency J2J is modified to 



Mnf 



1 Qinf ^crit ^eq j-j 

7) -"inf 

Z7T a cr it a e q a 



= 4.7 x 10 b GHz 



inf 



6 x 10- 5 M, 



pl 



3/4 



72 km/s • Mpc 



1/2 



0.1 mm 



1/4 



1 + Z, 



cq 



3200 



-1/4 



(22) 



Notice that the above estimate neglects another remarkable high-energy effect caused by the excitation of Kaluza- 
Klein modes (KK- modes). The KK-modes can propagate into the bulk and may be observed on the brane as massive 
gravitons. By contrast, the GW propagating on the brane is called the 'zero-mode', and it behaves like a massless 
graviton on the brane. Strictly speaking, the KK-modes and zero-mode are coordinate-dependent concepts and are 
mathematically well-defined only in the case of the Minkowski brane and the de-Sitter brane. Nevertheless, we keep 
to use these terms in the FRW case to distinguish these propagation features. 

As we will see in the next section, the brane is generally moving in the bulk spacetime. From the analogy of 
the moving mirror problem |4r|. even if only the zero-mode is generated in the inflationary epoch, the zero-mode is 
partially transferred to KK-modes with the arbitrary masses. If this is true, the amplitude of the zero-mode observed 
on the brane may be suppressed in comparison with the result neglecting the KK-modes, i.e., (|21|) . The envisaged 
spectrum involving the two effects is schematically shown as solid lines in Fig. [21 Unfortunately, we cannot fully 
estimate the KK-mode effects in an analytic way. Thus, in order to construct the IGWB spectrum including the 
KK-mode effects, we must directly solve the evolution equations of GWs in the five-dimensional cosmology. 

In the next section, we present the formalism to solve the evolution equations for GWs numerically. 



MD RD(low-energy) RD(high-energy) 



10" 6 

a io- 10 



10 -14 



fh ft 



f 4 '; 
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FIG. 2: A schematic diagram of the spectrum of IGWB. The four-dimensional prediction is shown in long-dashed lines. 
Considering the modification due to the non-standard cosmological expansion, the spectrum behaves as J7gw oc / 4//3 shown 
in short-dashed line, which may appear upon the detection limit of advanced LIGO (dot-dashed line). Moreover, KK-mode 
excitations modify the spectrum as solid lines. The main issue in this paper is to clarify whether the resultant spectrum (solid) 
exceeds the four-dimensional prediction (long-dashed). 



III. EVOLUTION EQUATION AND INITIAL CONDITIONS 
A. Evolution equation of GWs 

Let us consider the tensor perturbations in the AdSs spacetime. In the Poincare coordinate (r, x, z), the perturbed 
metric of the AdSs spacetime is given by 

ds 2 = f-\ {-dr 2 + (Sij + h ij )dx i dx j + dz 2 }, = 1, 2, 3), (23) 
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where hij satisfies the transverse-traceless (TT) condition, 

ditfj = h\ = 0. (24) 

While this coordinate system is free from the coordinate singularities, the brane is non-static, moving in the AdSs 
spacetime |4ll l42j|. The trajectory of the moving brane is determined from the scale factor on the brane, which is 
described as (rb, Zb): 

r h =T(t), Zh = ^)' ( 25 ^ 

where we use the cosmic time t on the brane as a parameter of the trajectory. The function T(t) is given by |47^ 

1 



T(t) = -y/T+iMf. (26) 



a 



In Fig. |2| we show the trajectory denoted by 'Friedmann brane' in the conformal chart where the surfaces of r = 
const, and z = const, are plotted as dashed lines. One can check that this trajectory induces the metric of the 
four-dimensional flat Friedmann- Robertson- Walker model on the brane : efeb = —dt 2 + a 2 (t)8ijdx t dx 3 . Note that, in 
the case of de Sitter brane, the scale factor becomes a(t) — e Ht and the Hubble parameter H is constant. Hence the 
trajectory becomes 



H 



= le~ n \ (27) 



which yields the straight trajectory in the bulk; that is, t£ s becomes proportional to z{J s . 

Let us focus on the evolution of the tensor perturbations hij. For convenience, we decompose the quantity hij into 
the three-dimensional spatial Fourier modes as 

Mr,x,*) =J2 I h^(r,z)e ik ^d 3 k, (28) 
p J 

where denotes a transverse-traceless polarization tensor. In terms of the Fourier modes, the evolution equation 
for the perturbations is given by |4S| 

, d 2 h d 2 h 3dh - „ 
O A dSsh=—- 1 — + -— + k 2 h = 0., 29 

ot z az z z oz 

where DAdSs denotes the d'Alembertian operator in the AdSs spacetime. Hereafter we omit subscripts of and 
simply write h. The equation l|29|l must be solved with the junction condition imposed on the brane. In the Poincare 
coordinate, the explicit form of the junction condition becomes |48j 



dh 
dn 



branc 



yi + H 2 e 2 dh 

Hi ~dz 



= 0, (30) 

z=z b (t) 



where d/dn denotes the normal vector of the brane. While there is generally a contribution from the tensor part of 
anisotropic stress tensor on the brane 7rS, we neglect it for simplicity and set ttJj = hereafter. 

It is important to note that the evolution equation (|29|l can be written in a separable form and by using this fact, 
one obtains the general solutions (e.g. [2J, \ 



h(T,z) = J dm\L 1 {m)z 2 H i 2 1 \mz)e l ^ T + h 2 {m)z 2 H ( 2 2) (mz)e- luJT j, (31) 

where lo 2 = m 2 + k 2 . The functions and denote respectively the Hankel functions of first and second 

kind, and hi(m) and /i2(w) represent arbitrary coefficients. The above expression implies that the GWs propagating 
in the bulk are described as a superposition of the zero mode (m = 0) and the KK-modes (m > 0). Solving the 
wave equation l|29() with the junction condition l|30|l is the equivalent task to determining the coefficients hi^im) 
that satisfy the junction condition |4^. In the very high-energy case, technical difficulties hinder efforts to calculate 
these analytically because of the significant contribution from the massive modes (ml S> 1). For this reason, we solve 
numerically the wave equation l|29|l with the junction condition i|3Ufl . 
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FIG. 3: The motion of the de Sitter and Friedmann brane in the AdSs bulk in the Poincare coordinate system. 



B. Initial conditions 



In order to correctly estimate the effects of the KK-modes, we must specify the initial conditions for the perturbed 
quantity h after the inflation. In this paper, we specifically consider a brane inflation model in which the exponential 
expansion takes place on the brane. According to Ref. the GN coordinate system (t,x,y) provides a useful 

spatial slicing in the inflationary epoch. With this coordinate, the perturbed metric of the AdSs spacetime becomes 

ds 2 = ~N 2 {y, t)dt 2 + A 2 (y, t){5 l3 + h ij )dx i dx i + dy 2 , (32) 

where A(t,y) = e Ht N(y) and N(t,y) = N(y) — cosh(y/£) — (1 + p/X) smh(y/t). Note that the perturbations hij 
satisfies the TT conditions <|24ll . In this coordinate, our brane is located at a fixed point y = 0. 

During de Sitter inflation, the solution of the evolution equation of GWs, □AdSs' 1 = [see Eq. 1)29(1]. can be 
obtained in a separable form, h(t, y) = u(y)(f>(t). Introducing the separation constant to which represents the mass of 
KK-modes, the mode functions u and (f> satisfy the following equations |40j : 

3H-— + [m + -t ) <pm = 0, (33) 



dt 2 dt \ a 

d 2 u m . N' du. 



dy 2 ' 4 iV^f + A^ = °< ^ 

where a prime denotes a derivative with respect to y. Picking up the normalizable modes from the solutions of the 
equation one notice that a large mass gap arises between the lightest KK-mode (m = 3H/2) and the zero-mode 
(to = 0) [4(| • From the point of view of the quantum theory, the large mass gap highly suppresses the excitation of KK- 
modes during inflation. Consequently, the zero-mode solution in the GN coordinates gives a dominant contribution 
to the metric fluctuation |4jj. Solving the equations l|3lfl) and l|3*4"|) with to = 0, the zero-mode is described explicitly 
as 

h{i 1 ,y) = C(-ki 1 f' 2 H^ 2 {~kr 1 ) 1 (35) 

where C is a normalization constant and r\ is the conformal time r] = —1/aH. 

To see how the solution l|35|l looks like in the Poincare coordinates, we rewrite the general solution (|31|l in the 
GN coordinates. In the de Sitter case, the coordinate transformation between the GN coordinates and the Poincare 
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coordinates is explicitly given as |2J, |25j 

z = — rjsmh(y/£), t — r\ cosh(y/£), (36) 
Substituting them into the general solution H31JI . we obtain 

h(T,z) = J dm {rj sinh(y/^)} 2 [h^m) {mr, sinh(y/f)) e -*« 

+ fca(m) i^ 2) (m77sinh(y/^)) e *«"J«*fr(»/A j _ ( 37 ^ 

Comparing (|35|1 with (|37[1 . we see that the zero- mode solution given in the inflationary epoch cannot be simply 
expressed in terms of the zero-mode solution in the Poincare coordinates, which indicates that a mixture of KK- 
modes is required to express the zero-mode solution in the inflationary epoch. Nevertheless, in the long-wavelength 
limit k — > 0, both the zero-mode solutions become constant over the time and the bulk space, and they coincide 
with each other. Since we are specifically concerned with the evolution of long-wavelength GWs after inflation, the 
constant mode, i.e., h = const, and dh/dr = 0, seems a natural and a physically plausible initial condition for our 
numerical calculation in the Poincare coordinate. Strictly speaking, however, this is valid only in the long-wavelength 
limit k — > 0. This point will be discussed in details in Sec IV Bl 



IV. NUMERICAL SIMULATION 



A. Numerical scheme 



On the basis of the formalism presented in the previous section, we now discuss the numerical treatment used to 
solve the wave equation l|29l) with the junction condition H30|) . First of all, the computational domain should be finite. 
We introduce an artificial cutoff (regulator) boundary in the bulk at z = z lcg (shown in Fig. [3J and impose the 
Neumann condition at the regulator boundary, i.e., 



dh 

dz 



0. (38) 



In this paper, numerical calculations were carried out by employing the pseudo-spectral method |49| . The amplitude 
of GWs h(r, z) is decomposed in terms of Tchebychev polynomials, defined as 

T„(£) ee cos^cos" 1 ^)) for - 1 < f < 1, (39) 
which yields a polynomial function of £ of order n. For example, 

2b(0 = l. r i(0=£, T 2 (£)=2£ 2 -l,... (40) 

Here the variable £ is related to the Poincare coordinate z. To implement the pseudo-spectral method, instead of 
using the Poincare coordinates (t,z) directly, we use the new coordinates (t, £): 

r = T(t), z=~[{z ieg -z b (t)}£ + {z Tes + z h {t)}] (41) 

so that the locations of both the physical and the regulator branes are kept fixed, and the spatial coordinate z 
is projected to the compact domain — 1 < £ < 1. Adopting this coordinate system, the amplitude h(t,£) is first 
transformed into the Tchebychev space through the relation, 

N 

h(t,0 = Y,hn(t)T n (0, (42) 

re=Q 

where we set N — 2048 or 4096. We then discretize the £-axis to the N + 1 points (collocation points) using the inho- 
mogeneous grid = cos(nir/N) called Gauss-Lobatto collocation points. With this grid, fast Fourier transformation 
can be applied to perform the transformation between the amplitude h(t, £) and the coefficients h n (t). Then the wave 
equation (|29|l rewritten in the new coordinates are decomposed into a set of ordinary differential equations (ODEs) 
for h n (t). For the temporal evolution of h n (t), we use Adams-Bashforth-Moulton method with the predictor-corrector 
scheme. Further technical details of the numerical scheme is summarized in Appendix lAl 
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B. Setup and parameters 

We are especially concerned with the late-time evolution of GWs after the inflation. For this purpose, we focus on 
the evolution equation l|29|) in the RD phase. In order to quantify high-energy effects, we define a useful parameter 
which represents the normalized energy density at the horizon re-entry time of the GWs concerned, namely, 
e* = e(t*)- From the Friedmann equation l|14|) and the definition of the scale factor <|16f) . the comoving wave number 
of GWs is rewritten in terms of the parameter e* as 



k = a*H* = (^-j v / e* + 2e *- (43) 

For higher frequency GWs, e* becomes larger. One thus expects that the high-frequency GW modes tend to be 
significantly affected by the high-energy effects. 

Notice that the location of the regulator brane is another important parameter. Here, the location of the boundary 
is set to z rog = 25-200£, which is far enough away from the physical brane to avoid artificial suppression of light 
KK-modes. Further, we must stop the numerical calculations before the influence of the boundary condition at 
z = z rcg reaches the physical brane z^. The arrival time of the influences of the regulator brane can be estimated by 
drawing a null line from the initial position of the regulator boundary (tq, z rcg ) toward the physical brane. With these 
treatments, we have checked that the amplitude of GWs on the brane is fairly insensitive to the location of regulator 
boundaries. Thus, all the results presented in Sec. Ivland Sec. I VII are free from the effect of regulator boundary. 

In the situation considered here, the initial time i; n j t is also an important parameter, which turns out to have an 
important effect on the GWs in the bulk |3(j • We parameterize the initial time as 

a(tinit)H(t init ) 

•Sinit — ; V"7 

which represents the wavelength of GWs normalized by the Hubble horizon scale at the initial time i; n it- In order to 
get a reliable estimate, we set Sj n i t 3> 1 and run the simulations until e(t) <^ 1, when the high-energy effects on the 
GWs become negligible. 

Finally, we adopt the constant mode /i(ii n it,£) = 1 as an initial condition according to the discussion in Sec. IIII Bl 
Although it is plausible, the validity of the constancy of the superhorizon modes must be checked. This point will be 
carefully discussed in Sec. IV Bl 



l/3(l+u>) 



V. CODE CHECK AND QUALITATIVE BEHAVIOR OF GWS 



A. Code check 



In order to check our numerical code, we first consider the simplest case in which the location of the brane is given 
by z — Zb =constant. This is the so-called Minkowski brane embedded in the AdS5 bulk. The solution of the evolution 
equation is given as plIM!^ 

fr ex act(T,z) = z 2 Z 2 (mz)E(ujT) , (45) 

where u> = Vk 2 + m 2 . The function Z 2 and E respectively denote linear combinations of the Bessel functions of order 
2 and the sinusoidal functions. Imposing the junction condition (|30fl at Zb = we obtain 

ftexact(T,z) = {Y 1 (m)J 2 (mz) - J 1 (m)Y 2 {mz)} cos(ivt) . (46) 

According to the analytic solution (|46|) . we set h = /i OX act(0, z) and h = <9/i GX act / 9t\ t= q = for the initial condition 
of numerical simulation and compare the numerical results with the analytic solution. 

Fig. 21 shows the behavior of GWs in the AdSs bulk. The right panel of the figure is the projection of the left 
panel. In this simulation, we chose parameters as z rcg = 20^?, ml = kl = 2 and N = 1024. The left panel of Fig. [S] 
shows the snapshot of the waveform at r = T\ = 10£, which illustrates that the numerical result accurately recovers 
the exact solution (|46|) in the interval between Zb < z < z rcg — t\ . Outside this region, the numerical simulation is 
contaminated by the boundary condition of the regulator brane. In the right panel of Fig. [S] the fractional error of 
the amplitude \(h num — h cxa _ ct ) / h oxact \ evaluated at the time r = t\ is plotted as a function of bulk coordinate. We 
found that the error is suppressed to the order of 10~ 3 near the physical brane. Note that there appear several sharp 
spikes, whose locations roughly match those of the zero-point ft.oxact(i", z) — 0. Thus, the cancellation of significant 
digits occurs. These numerical errors can be reduced when the number of collocation points N is increased. 
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FIG. 4: The behavior of a test wave with ml = k£ = 2 in the bulk. The Minkowski brane is located at z = I. The right panel 
depicts the projection of the three-dimensional waves of the left panel. 
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FIG. 5: The left panel illustrates that the numerical solution at r = 10 is consistent with the analytical solution l|46|l in the 
region Zh < z < 10. The outer region 10 < z < z Icg in the bulk is contaminated by the boundary condition on the regulator 
brane at z = z ICS - The right panel shows the numerical errors |(/i n um — /icxact)//iexact| estimated at that time, which is suppressed 
by 10~ 3 near the physical brane. The spike-shapes reflects the cancellation of significant digits because of /i eX act(10, z) w 0. 



B. Behavior of GWs in the bulk and the validity check of the initial condition 

Having checked the reliability of our numerical scheme, we now focus on the cosmological evolution of GWs. As 
we mentioned in Sec IIII Bl we must first clarify the validity and the sensitivity of the initial condition, namely, the 
constancy of the superhorizon modes. It should be stressed that, only in the long-wavelength limit k — > during the 
inflationary phase, the constant mode coincides with the zero- mode solution [see Ea.(|37|)], Therefore, the constancy 
of GW amplitudes after inflation cannot be guaranteed even on superhorizon scales. Depending on the choice of the 
parameters Si n it and e*, the mode h = const, may not be a good approximation to the initial condition for numerical 
simulations in the RD epoch. 

Fig. shows the time evolution of GWs in the Poincare coordinate system in the de Sitter case. In this simulation, 
we set Si n it = 100 and HI = V3. The universe on the brane experiences accelerated cosmological expansion and the 
wavelength of GWs becomes longer than the Hubble horizon. Fig. [(^indicates that the GWs on initially superhorizon- 
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scale remains constant not only on the brane but also in the bulk. The right panel of this figure, which depicts the 
projection of the left panel, shows that a very slight change of the amplitude is observed (a fraction of the original 
amplitude of < 1%) and the amplitude finally converged to a fixed value. In this sense, the constant mode h = const. 
is suitable for the initial condition of the superhorizon-scale GWs in the RD phase if we impose the initial condition 
just after the inflation. 

Adopting this initial condition, we then performed simulations in the radiation-dominated FRW case (w = 1/3) 
with same parameters, e* = 1.0 and Sj n i t = 100. This result is shown in Fig d We found that the constancy of the 
GW amplitude no longer holds in the bulk even before r *=s 0, where the wavelength of GW on brane just re-enters 
the Hubble horizon. In particular, GWs emanating from the physical brane are observed, which propagate into the 
bulk spacetime almost along a null line. This indicates that the excitation of KK-modes occurs near the brane even 
if the wavelength of GWs is still outside the Hubble horizon. 




zll 



FIG. 6: The evolution of a GW in the bulk in the case of a de Sitter brane. We set the Hubble parameter to Hi — v3 with 
(sinit, -z r og) — (10,20). The right panel depicts the projection of the three-dimensional waves of the left panel, zooming in the 
image in < z/i < 3. The empty corner in the surface represents the motion of the brane [see Ea. l|27|l ], 




FIG. 7: The evolution of a GW in the bulk in the case of a Friedmann brane. We set the comoving wave number to k — y/3/i 
or e* = 1.0 with (smit, Zrcg) = (200, 80). The right panel depicts the projection of the three-dimensional waves of the left panel. 
The empty corner in the surface represents the motion of the brane [see Eq. l25H l. 

It is noteworthy that the different behaviors of GWs in the AdSs bulk may be caused by the difference in the motion 
of the brane [see Eas. i|25[l and i|27|) ]. In the moving mirror problem in an electromagnetic field, the acceleration or 
deceleration of the mirror yields the creation of photons due to vacuum polarization (e.g., Sec. 4.4 of Ref. p^h 
A similar phenomenon may occur in the AdSs bulk, that is, the KK-modes (massive gravitons) are excited by the 
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FIG. 8: Snapshots of the GW amplitudes in the bulk for various choices of initial time. The snapshots were taken when the 
wavelength of GWs becomes five times longer than the Hubble horizon, i.e., aH/k = 5. 
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FIG. 9: Evolved results of GWs projected on the brane starting with the various initial times. 



deceleration of the brane which is depicted by an arc with a non-zero curvature cPzb/dr 2 < 0. 

Figs. El and [7| reveal that the constant mode /i(ii n it,£) = const, can be used as the initial condition if we set this 
just after the end of inflation, but, the constancy of the long-wavelength mode would not be guaranteed in the RD 
epoch even on super horizon scales. This implies that the choice of the initial time ti n n (or s; n it) defined in (|44ll is 
crucial when setting the initial condition at the RD epoch. 

In Figs. OH and EH the dependence of the evolution of GWs on the initial time is shown by varying the parameter 
Si n it in low-energy (e* = 0.1) and high-energy (e* = 10) cases. Fig. plots the snapshots of the amplitude h(r 7 z) in 
the bulk when the wavelength of GWs becomes five times larger than the Hubble horizon, i.e., aH/k = 5. Clearly, in 
the bulk, the amplitude of GWs is very sensitive to the choice of the parameter Sj n j t , or equivalently, the initial time 
tinit- The resultant wave- form away from the physical brane does not show any convergence even in the low-energy 
case (e* = 0.1). This behavior may be caused by the fact that the constant mode with the comoving wave number k 
in the RD epoch immediately starts to oscillate as h(r, z) oc e lkT , which is the massless (m — ► 0) limit of Eq. I|31l) . 

On the other hand, in Fig. [51 the GW amplitudes tend to converge on the brane if we set the initial time i; n ;t early 
enough. This convergence property might be due to the presence of the junction condition (|30() . Therefore, as far as 
we choose Si n it ^ 50 for our interest of the energy scale 0.01 < e* < 100, we do not need to care about the initial 
time, when we estimate the IGWB spectra on the brane. In Appendix [51 quantitative aspects of the convergence 
properties of the amplitude are discussed. Moreover, Seahra addressed these points in an analytic way in Ref . |34| . 
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FIG. 10: Squared amplitude of GWs on the brane in low-energy (left) and the high-energy (right) regimes. In both panels, 
solid lines represent the numerical solutions of wave equation 12911 . The dashed lines are the amplitudes of reference wave ft re f 
obtained from equation 14711 . 



VI. IGWB SPECTRA 



A. Comparison with reference waves 

Keeping the results in Sec. E]in mind, let us now quantitatively estimate the high-energy effects of the GWs and 
evaluate the energy spectra of the IGWB on the brane. To quantify these, it might be helpful to discriminate the 
influence of KK-mode excitation in the bulk from the non-standard cosmological expansion caused by the p 2 -term. 
in the Friedmann equation (|14|) . For this purpose, we introduce the reference wave h le f, which is a solution of the 
four-dimensional evolution equation of the amplitude obtained by replacing the scale factor and the Hubble parameter 
derived from the standard Friedmann equation with those from the modified Friedmann equation (|14f) . The resultant 
equation is given by 

fc^ 2 



h Te f + 3Hh Te f + href = 0, (47) 

which is just the Klein-Goldon equation for a scalar field in the FRW background (e.g. 0,0) and is same as for 
m = 0. Comparing the numerical simulations with the solution of the wave equation (|47|l . the effect of the KK-mode 
excitation can be quantified separately. 

Fig. EH shows the squared amplitude of the GWs, ftJ D and /i 2 ef as functions of the scale factor a. The left panel 
shows the low-energy case (e* = 0.1), while the right panel depicts the result in the high-energy regime (e* = 50). 
As we increase the energy scale at the horizon re-entry time, the GW amplitude h^u becomes considerably reduced 
compared to the reference wave, h ra f. Since the late-time evolution of GWs simply scales as h oc 1/a in both /15D and 
h rc [ , the suppression of the amplitude /15D is caused by the excitation of KK-modes around the horizon re-entry time. 
Notice that the normalized energy density at the horizon re-entry time e* is related to the observed proper frequency 
2ir f = k/a = (a*/a )H* as 

7^= — W*= — (48) 

Jcrit \ a crit/ \ £* / 

where the critical frequency / cr it is defined in (|20|l as 27r/ cr i t = (a cr it/ao)^ _1 , and w — 1/3 in this case. From this 
relation, one expects that the KK-mode excitation is essential in the high-energy regime and the deviation from the 
standard four-dimensional prediction for the spectrum of IGWB becomes more prominent above the critical frequency, 

/ > /crit- 
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FIG. 11: Frequency dependence of the ratio of amplitudes |/i5d//w | between the numerical simulation of wave equation 12911 
and the reference wave 1471 . The vertical solid line represents the critical frequency. The dashed line indicates the fitting result 
1491 . where fitting was examined using the data with Sinit = 200 at the asymptotic region e* > 5. 
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TABLE I: Numerical parameters used for the simulations to estimate the frequency dependence of the ratio of amplitudes 
\h5T>/h Te {\. 



B. IGWB spectrum in the five-dimensional cosmology 



We are in position to estimate the influence of KK-mode excitation on the shape of the spectrum. To do so, we 
ran simulations for the parameters listed in Table |Tj(io = 1/3 case) and estimated the ratio of amplitudes \h^/h le {\ 
for a different set of parameters. Note that in the simulations with e* < 1, the location of regulator brane z rC g should 
be set far away from the physical brane. This is because the long-term evolution is needed to follow the oscillatory 
behavior. 

We show the frequency dependence of the ratio in Fig. ^] The ratio is evaluated at the low-energy regime long 
after the horizon re-entry time and is plotted as a function of the frequency /// cr it- Clearly, the ratio \hsD/h Te f\ 
monotonically decreases with the frequency and the suppression of amplitude becomes significant above the 
critical frequency / cr it- Using the data points in the asymptotic region e* > 5, we try to fit the ratio of amplitudes 
with Si n it = 200 to a power-law function. Employing the least-squares method, the result becomes 



*5D 



J_ 

/crit 



(49) 



with a = 0.76 ± 0.01 and j3 = 0.67 ± 0.01 (dashed line in Fig. Ill|) . In Appendix iBl we calculate the ratios for various 
Sinit for each combination (e*, z reg ) to check the robustness of this result. 

The power-law fit (|49(l can be immediately translated to the energy spectrum of IGWB, Oqw- The spectrum taking 
account of the KK-mode excitations is calculated as 



^GW — 
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(50) 
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FIG. 12: The energy spectrum of the IGWB around the critical frequency. The filled circles represent the spectrum caused by 
the non-standard cosmological expansion of the universe. Taking account of the KK-mode excitations, the spectrum becomes 
the one plotted as filled squares. Particularly, in the asymptotic region depicted in the solid line, the frequency dependence 
becomes almost same as the one predicted in the four-dimensional theory (long-dashed line). 



where we used the fact ficw ^ h 2 f 2 . As discussed in Sec. Ill B 21 if we neglect the effect of the KK-mode excitation, 
the spectrum becomes f2 re f oc / 4 / 3 [See Eq|^. Then, combining it with the result l|49|l . the IGWB spectrum becomes 
nearly flat above the critical frequency : 

^gw oc f°, (51) 

which is shown in filled squares in Fig. 1121 In this figure, the spectrum calculated from the results of the reference 
waves f2 rc f is also shown in filled circles. Note that the normalization factor of the spectrum is determined so as to 
be f^GW = 10~ 14 according to the constraint from the CMB observation. The short-dashed line and the solid line 
represent each asymptotic behavior in the high-frequency region. The spectrum taking account of the two high-energy 
effects seems almost indistinguishable from the standard four-dimensional prediction shown in long-dashed line in the 
figure. In other words, while the effect due to the non-standard cosmological expansion lifts up the spectrum, the 
KK-mode effect reduces the GW amplitude, which results in the same spectrum as the one predicted in the four- 
dimensional theory. Additionally notice that the amplitude taking account of the two effects near / « is slightly 
decreasing, which agrees with the results in our previous study for e* < 0.3 using the GN coordinates [2£|. 

At this point, however, it is unclear whether the results obtained here is generic or accidental for certain range of 
the model parameters. To clarify the cosmological dependence of the KK-mode excitations in more quantitative way, 
we next study the cases with different EOS parameter w. 



C. Dependence on equation of state 

To quantify the EOS dependence, we ran simulations for the MD case w — and the somewhat stiff matter case 
w = 1, which might be realized by introducing the kinetically driven scalar field (e.g. the quintessential inflation 
|5lL l52|). Varying the EOS parameter changes the acceleration of the brane. One naively expects that the different 
motion of the brane may suppress or enhance the KK-mode excitation. 

With the same procedure as in the previous subsection, we calculated the ratio of amplitudes \h 5 D/h ro [ | for various 
e* in the case of w = and w = 1. The results are summarized in Fig. 1131 where the horizontal axis represents 
\/l + H 2 £ 2 . Fitting the power-law function to all cases shown in the figure, we found that the ratios universally scale 
as 

= 5(1 + H 2 J 2 y°- 2A « 5(1 + H 2 J 2 )~ l /\ (52) 
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FIG. 13: EOS dependence on the amplitude ratio. The filled circles, squares and triangles show the cases with w = 0, 1/3 
and 1, respectively. The solid line, short-dashed line and long-dashed line show each fitting result above the critical frequency 
corresponding to (1 + H^l 2 ) 1 ^ 2 — \/2 depicted in the vertical dot-dashed line. 



where a = 0.75, 0.81 and 0.83 for w = 0, 1/3 and 1, respectively, which indicates that the quantity 5 may be related 
to w as a — — 0.155iu 2 + 0.241u; + 0.748 by simply fitting the quadratic function. An important point to emphasize 
is that this scaling property does not depend on the parameter w or the acceleration of the brane. Combining the 
scaling relation with (|50() , the IGWB spectra can be estimated as 

n GW = a 2 (l + if^2)-i/2 firef _ (53) 

Especially, in the high-frequency region / ^> / cr it> the prefactor of the right-hand side behaves as 

(l + Hlt 2 )- 1 ' 2 kHJk (54) 
from equations Q,© an d (|18[) . Combining the result (|19|) . the energy spectrum of the IGWB behaves as 

n GW cx f*£* for / > / crlt . (55) 
Owing to these calculations (|19|l and JSSJl, for w — (MD), we obtain 



V2 for / » /, 



Ugw { J n J JCrit ' (56) 

\f- 2 for /«/crit, 

and for w = 1 case, 

J/ 2 /5 for />/ crit , 
^gw oc <^ (57) 
[f L for / < / crit . 

These results imply that the spectrum generally changes from the four-dimensional prediction. Indeed, transforming 
the numerical results 1)52(1 to the energy spectra in the cases with w — and w = 1, the frequency dependence 
of the spectra including the two high-energy effects (solid lines) clearly differ from each four-dimensional prediction 
(long-dashed lines). 

Taking these results into account, one may conclude that the cancellation of the high-energy effects in the RD epoch 
is accidental and the KK-mode excitation dominates over the non-standard cosmological expansion when w > 1/3. 
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FIG. 14: The energy spectrum of the IGWB in the background EOS with w = O(left), and with w = l(right). The amplitude 
is normalized by the value of the four-dimensional prediction at the critical frequency. 



VII. CONCLUSION 



We have investigated the power spectrum of the IGWB in the five-dimensional cosmology based on the Randall- 
Sundrum model. In the braneworld scenario, the two high-energy effects affect the shape of the spectrum above the 
critical frequency / cr i t defined in i|2(J|) . One is the non-standard cosmological expansion on the brane caused by the 
high-energy correction of the Friedmann equation. The analytical estimate taking account of this effect reveals that 
the effect makes the spectrum steeply blue [see Eq. H21fl ]. By contrast, another important effect is the excitations of 
KK-modes which escapes from our brane into the five-dimensional bulk, leading to the suppression of the spectrum. 
In order to quantify these two effects, we solved the wave equation of each Fourier mode of GWs numerically for 
various EOS parameters w. 

The systematic survey of numerical simulations with various parameter sets reveals that there may exist the universal 
scaling law for the KK-mode excitation in the high-energy regime [Ea. i|52[l ]: 

cx {i + H 2 j 2 y l /\ 

Using the universal scaling law, we constructed the power spectrum of the IGWB in the cases with w — (MD 
universe), w = 1/3 (RD universe) and w = 1 (stiff matter dominant universe). From the results l|12|) and (|55l) . the 
frequency dependence of the spectrum in the high-frequency region / > / cr ;t becomes 

{f- 2 (4D) , J- 1 ' 2 (5D) for w = 0, 
f°, f° for w = 1/3, (58) 

f 1 , J 2 / 5 for w = l. 

Particularly, in the RD case, the accidental cancellation of the two high-energy effects occurs, which yields the same 
spectrum as one predicted in the four-dimensional (4D) theory. This scaling law might be understood in a context of 
moving mirror problems. The discussion about the analytic derivation of the scaling law is work in progress. 

Finally, we briefly comment on the other numerical works using the different schemes. Recently, Seahra solved 
numerically the wave equation in the null coordinate system based on the Poincare coordinates using a sophisticated 
numerical scheme |34| . He observed the agreement between his numerical results and our results for the same choice 
of the EOS parameters with the same initial conditions as ours. In addition, it is observed that, even if a constant 
initial condition on the initial null hypersurface is chosen, the amplitude of GWs on the brane is almost identical with 
our results. Moreover, the numerical calculation based on the quantum theory has been performed by Kobayashi and 
Tanaka |33|| . They reported the same spectrum in the RD case as ours, even if KK modes are taken into account in 
the initial de Sitter phase. On the other hand, Ichiki and Nakamura have obtained a tilted spectrum flow oc f~ 0A6 
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[3l| . While their early results have included errors associated with numerical accuracy, the new calculation using the 
revised code did not converge to the flat spectrum either. Currently, we do not know the reason why the result by 
Ichiki and Nakamura is different from ours. In order to understand these numerical results well, the analytical study 
of the scaling relation 1|52|) is essential, which is definitely our next task. 
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APPENDIX A: THE PSEUDO-SPECTRAL METHOD 

In this appendix, we briefly describe the implimentation of the pseudo-spectral method in our numerical scheme. 
From the coordinate transformation l|41|l. the wave equation (|29|l is expressed as 

d 2 h d 2 h d 2 h dh dh , k N 

W + K «m + K «W + Kt ai + K6 di + Kh = ' (A1) 

where the coefficients K t ^, K^, K t , and K are functions of t and £. We use the predictor-corrector method for the 
temporal evolution. To implement this, we introduce an auxiliary variable x(i, £) satisfying the equation 

dh dh 

^=X-K t ^=F( X ,h';t^ n ), (A2) 

where the prime denotes the derivative with respect to £ [see En.([41jl]. With this definition, the time evolution of \ 
satisfying to H29[) is formally written as 

^ = G( X ,h,h',h";t^ n ). (A3) 
at 

Notice that the function G docs not contain the derivative <9x/<9£. Empirically, the presence of this derivative causes 
numerical instability. The functions F and G are evaluated at each collocation point £„. Then, transforming to the 
Tchebychev space by Eq. (|42[1 . we obtain a set of ordinary differential equations : 

ir = p »®> % = dn{t) - (A4) 

With the reduced equations l|A4|l . the predictor-corrector method based on the Adams-Bashforth-Moulton scheme 
can be used to obtain the time evolution of h n (t) and x~ n (t). 

At each time step, while the boundary conditions (|3Ufl and (|38|l are used to evaluate h^-i and hpj, we put additional 
conditions on Xn-i an d Xn as 

Xn-i = Xn = 0. (A5) 

Owing to the definition of the function F containing no derivatives of x [ see Eq. (|A2ll ]. this empirically based treatment 
of boundary conditions suppresses the numerical error caused by finite truncation of the Tchebychev transformation 

In summary, we firstly evaluate the functions F and G in the physical space at the time t. Then, transforming 
them into the Tchebychev space by l|42(l . we obtain h n and Xn at the next time step t + At by solving the reduced 
equations 1A4(I for < n < N — 2, and imposing the boundary conditions and the additional conditions l|A5() into h n 
and Xn for n = N — 1, AT. The spectral coefficients of the derivatives h' and h" can be computed in decreasing order 
by the recurrence relations 

c J« = h% + 2(n+ l)h n+1 , h% N = 0, (A6) 

and then 

c~h^ = h% + 2(n+ 1)/£> X , h% N - 0, (A7) 
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FIG. 15: Convergence of the ratio of amplitudes /isd/Zw- Solid lines represent the fitting curve calculated by employing the 
nonlinear least-square method. Long-dashed lines denote the asymptotic value of the ratio. 



where c„ is defined as 



2 for n = 0, TV 

1 for 1 < n < N - 1 



(A8) 



Finally, we obtain h(t + At, £) and its derivatives h'(t + At, £) and h"(t + At, £) by the inverse Tchebychev transfor- 
mation. 



APPENDIX B: INITIAL TIME DEPENDENCE 



As seen in IV 51 the amplitudes on the brane becomes insensitive to the choice of the initial time as long as Smit is 
large enough. In this appendix, we show that the ratios of amplitudes discussed in IVI Bl tend to converge to a fixed 
value in the low-energy and high-energy cases. This validates the estimation of the power spectrum of the IGWB 
using the results obtained from the simulations with Si„;t = 200. 

Fig. 1151 shows the dependence of the ratios /iso/^ref on the parameter s; n it. The left and right panels show the 
high-energy (e* = 50) and the low-energy (e* = 0.1) cases, respectively. In both cases, the ratios clearly converge to 
certain asymptotic values as increasing Si n it- Using the nonlinear least-square method, we tried to fit these values to 
the function 



*5D 



h Te 



— Sib; 



a 



where A,B,C are fitting parameters depending on e*. Then we obtained 



(A,B,C) 



(0.566, -0.573, 0.772) for e* = 0.1, 
(0.456, -0.880, 0.118) for e* = 50, 



(Bl) 



(B2) 



which are shown in solid curves in Fig. (|15|) . Note that C represents the asymptotic values shown in long-dashed 
lines in each panel. 

Picking up the values at s; n it = 200 in both cases, one can see that the deviations from the asymptotic values C 
keep to be less than a few percent. From this fact, we use the ratios with Si n ; t = 200 to construct the power spectra 
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of the IGWB without deriving the asymptotic values for each e* . 
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